module valuefunction

use params
use states
use searchnew
use wage
use utilityfs

implicit none

contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine solvevalfun_deadparent(nodek)

  implicit none

  real(prec) cashkNT, cashkPT, cashkFT
  real(prec) valNT, valPT, valFT
  real(prec) leisureNT, leisurePT, leisureFT
  real(prec) tempconNT, tempconPT, tempconFT
  real(prec) ax, bx, cx

  integer nodek
  integer wknode,aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  ax=0.0d0

  !1. not working
  leisureNT=get_leisure(workhrs(1),famhrs(1))
  cashkNT=max(assetk(aknode,age)+0.0*incomek(wknode,age)+nonwageinc(wknode),cmink)
  cx = cashkNT ! max possible expenditure (no borrowing)
  bx = cx*0.6d0
  valNT = golden_cons_deadparent(ax,bx,cx,valfun_k_deadparent,tempconNT,leisureNT,cashkNT,nodek)
  valuekNT_deadparent(age,nodek) = -valNT !in golden we solve for min of negative objective func
  conexpkNT_deadparent(age,nodek) = tempconNT

  !2. part-time
  leisurePT=get_leisure(workhrs(2),famhrs(1))
  cashkPT=max(assetk(aknode,age)+0.5*incomek(wknode,age)+nonwageinc(wknode),cmink)
  cx = cashkPT ! max possible expenditure (no borrowing)
  bx = cx*0.6d0
  valPT = golden_cons_deadparent(ax,bx,cx,valfun_k_deadparent,tempconPT,leisurePT,cashkPT,nodek)
  valuekPT_deadparent(age,nodek) = -valPT !in golden we solve for min of negative objective func
  conexpkPT_deadparent(age,nodek) = tempconPT

  !3. full-time
  leisureFT=get_leisure(workhrs(3),famhrs(1))
  cashkFT=max(assetk(aknode,age)+1.0*incomek(wknode,age)+nonwageinc(wknode),cmink)
  cx = cashkFT ! max possible expenditure (no borrowing)
  bx = cx*0.6d0
  valFT = golden_cons_deadparent(ax,bx,cx,valfun_k_deadparent,tempconFT,leisureFT,cashkFT,nodek)
  valuekFT_deadparent(age,nodek) = -valFT !in golden we solve for min of negative objective func
  conexpkFT_deadparent(age,nodek) = tempconFT

  !choose between labor supply options
  if (valFT .lt. valNT .and. valFT .lt. valPT) then
    LFP_deadparent(age,nodek) = 3
    valuek_deadparent(age,nodek) = valuekNT_deadparent(age,nodek)
    conexpk_deadparent(age,nodek) = tempconFT
  else if (valPT .lt. valNT .and. valPT .lt. valFT) then
    LFP_deadparent(age,nodek) = 2
    valuek_deadparent(age,nodek) = valuekPT_deadparent(age,nodek)
    conexpk_deadparent(age,nodek) = tempconPT
  else
    LFP_deadparent(age,nodek) = 1
    valuek_deadparent(age,nodek) = valuekFT_deadparent(age,nodek)
    conexpk_deadparent(age,nodek) = tempconNT
  endif

end subroutine solvevalfun_deadparent
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
real(prec) function valfun_k_deadparent(conexp,leisure,cashk,nodek)

  real(prec) conexp
  real(prec) leisure
  real(prec) cashk
  real(prec) saving
  real(prec) utnow
  real(prec) valtom, valtombit
  real(prec) cashtom, wgt1,wgt2,node_interp,estim1,estim2
  real(prec) valterm
  integer mina, nodeto, ww  

  integer nodek, wknode, aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  wgt1=0.0
  wgt2=0.0

  utnow = utility_deadparent(conexp,leisure)

  if (ageto>=nt+1) then
    valterm=utilitykterminal(cashk-conexp,wknode)
	valfun_k_deadparent = -(utnow+valterm)
  else

    !if not in final period, find next period asset stock on grid & calculate expected value (given alive next period)
    saving = r*(cashk - conexp)
    if (saving.le.minass) saving = minass

    call find_ass(2, saving, ageto, mina, wgt1,wgt2)
    node_interp=wgt1*real(mina)+wgt2*real(mina+1)

    !loop through future income states
    valtom = 0.d0
    do ww = maxwn(wknode,ageto)+1 ,minwn(wknode,ageto)-1

      !define income*asset (compact) node
      nodeto = (ww-1)*na

      !linear interpolation (not tri-linear cuz only interpolate on own assets)
      valtombit=interp_deadparent(ageto,nodeto,mina,node_interp,valuek_deadparent) !1st=insurance,3rd=health (both placeholders, dont matter for kid)
      valtom = valtom + valtombit*wprobk(ww,wknode,ageto)   !Weight by prob of ww at t+1 given wnode at t

    enddo   !end loop through future income states

    valfun_k_deadparent = - (utnow + (disc**(-1.d0))*(valtom))

  endif

end function valfun_k_deadparent
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!solve value function in divorce
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine solvevalfun_divorced(nodek) !no option of family care here
  
  implicit none
  
  real(prec) valk_div,valp_div
  real(prec) conk_div,conp_div
  integer LFP_div
  real(prec) maximum_valk_temp, maximum_conk_temp, maximum_asschoicek_temp
  real(prec) maximum_valk, maximum_conk, maximum_asschoicek
  real(prec) maximum_valp, maximum_conp, maximum_asschoicep

  integer i,acp,ack ! counters

  integer nodek !kid asset node

  !run through insurance / no insurance
  do i=1,ni

    LFP_div=1
    maximum_valp=-99999999.0

    do acp=1,nachoice !loop through parent asset choice
      
      maximum_valk_temp=-99999999.0
      
      do ack=1,nachoice !loop through kid asset choice

        call solvevalfunk_divorced(achoice(acp),achoice(ack),i,valk_div,LFP_div,conk_div,nodek)

        !kid chooses max kid asset for each parent asset
        if (maximum_valk_temp<valk_div) then
          maximum_valk_temp=valk_div
          maximum_conk_temp=conk_div
          maximum_asschoicek_temp=achoice(ack)
        endif

      end do


      !solve parent problem: 
      call solvevalfunp_divorced(achoice(acp),maximum_asschoicek_temp,maximum_conk_temp,LFP_div,i,valp_div,conp_div,nodek) !depends on kid value

      !parent chooses max parent asset
      if (maximum_valp<valp_div) then
          maximum_valp=valp_div
          maximum_conp=conp_div
          maximum_asschoicep=achoice(acp)
          maximum_valk=maximum_valk_temp
          maximum_conk=maximum_conk_temp
          maximum_asschoicek=maximum_asschoicek_temp
      endif

    end do

    valuepNT_divorced(age,i,bignode,apnode,nodek)=maximum_valp
    valuekNT_divorced(age,i,bignode,apnode,nodek)=maximum_valk
    conexppNT_divorced(age,i,bignode,apnode,nodek)=maximum_conp
    conexpkNT_divorced(age,i,bignode,apnode,nodek)=maximum_conk

    !***************************************************

    !***************************************************
    LFP_div=2
    maximum_valp=-99999999.0

    do acp=1,nachoice !loop through parent asset choice
      
      maximum_valk_temp=-99999999.0
      
      do ack=1,nachoice !loop through kid asset choice

        call solvevalfunk_divorced(achoice(acp),achoice(ack),i,valk_div,LFP_div,conk_div,nodek)

        !kid chooses max kid asset for each parent asset
        if (maximum_valk_temp<valk_div) then
          maximum_valk_temp=valk_div
          maximum_conk_temp=conk_div
          maximum_asschoicek_temp=achoice(ack)
        endif

      end do

      !solve parent problem: 
      call solvevalfunp_divorced(achoice(acp),maximum_asschoicek_temp,maximum_conk_temp,LFP_div,i,valp_div,conp_div,nodek) !depends on kid value

      !parent chooses max parent asset
      if (maximum_valp<valp_div) then
          maximum_valp=valp_div
          maximum_conp=conp_div
          maximum_asschoicep=achoice(acp)
          maximum_valk=maximum_valk_temp
          maximum_conk=maximum_conk_temp
          maximum_asschoicek=maximum_asschoicek_temp
      endif

    end do

    valuepPT_divorced(age,i,bignode,apnode,nodek)=maximum_valp
    valuekPT_divorced(age,i,bignode,apnode,nodek)=maximum_valk
    conexppPT_divorced(age,i,bignode,apnode,nodek)=maximum_conp
    conexpkPT_divorced(age,i,bignode,apnode,nodek)=maximum_conk
    !***************************************************

!     !***************************************************
    LFP_div=3
    maximum_valp=-99999999.0

    do acp=1,nachoice !loop through parent asset choice
      
      maximum_valk_temp=-99999999.0
      
      do ack=1,nachoice !loop through kid asset choice

        call solvevalfunk_divorced(achoice(acp),achoice(ack),i,valk_div,LFP_div,conk_div,nodek)

        !kid chooses max kid asset for each parent asset
        if (maximum_valk_temp<valk_div) then
          maximum_valk_temp=valk_div
          maximum_conk_temp=conk_div
          maximum_asschoicek_temp=achoice(ack)
        endif

      end do

      !solve parent problem: 
      call solvevalfunp_divorced(achoice(acp),maximum_asschoicek_temp,maximum_conk_temp,LFP_div,i,valp_div,conp_div,nodek) !depends on kid value

      !parent chooses max parent asset
      if (maximum_valp<valp_div) then
          maximum_valp=valp_div
          maximum_conp=conp_div
          maximum_asschoicep=achoice(acp)
          maximum_valk=maximum_valk_temp
          maximum_conk=maximum_conk_temp
          maximum_asschoicek=maximum_asschoicek_temp
      endif

    end do

    valuepFT_divorced(age,i,bignode,apnode,nodek)=maximum_valp
    valuekFT_divorced(age,i,bignode,apnode,nodek)=maximum_valk
    conexppFT_divorced(age,i,bignode,apnode,nodek)=maximum_conp
    conexpkFT_divorced(age,i,bignode,apnode,nodek)=maximum_conk
    !***************************************************

    !choose labor supply
    if (valuekFT_divorced(age,i,bignode,apnode,nodek) >= valuekPT_divorced(age,i,bignode,apnode,nodek) .and. &
        valuekFT_divorced(age,i,bignode,apnode,nodek) >= valuekNT_divorced(age,i,bignode,apnode,nodek)) then
      valuek_divorced(age,i,bignode,apnode,nodek) = valuekFT_divorced(age,i,bignode,apnode,nodek)
      valuep_divorced(age,i,bignode,apnode,nodek) = valuepFT_divorced(age,i,bignode,apnode,nodek)
      conexpk_divorced(age,i,bignode,apnode,nodek) = conexpkFT_divorced(age,i,bignode,apnode,nodek)
      conexpp_divorced(age,i,bignode,apnode,nodek) = conexppFT_divorced(age,i,bignode,apnode,nodek)
      LFP_divorced(age,i,bignode,apnode,nodek) = 3
    else if (valuekPT_divorced(age,i,bignode,apnode,nodek) >= valuekFT_divorced(age,i,bignode,apnode,nodek) .and. &
             valuekPT_divorced(age,i,bignode,apnode,nodek) >= valuekNT_divorced(age,i,bignode,apnode,nodek)) then
      valuek_divorced(age,i,bignode,apnode,nodek) = valuekPT_divorced(age,i,bignode,apnode,nodek)
      valuep_divorced(age,i,bignode,apnode,nodek) = valuepPT_divorced(age,i,bignode,apnode,nodek)
      conexpk_divorced(age,i,bignode,apnode,nodek) = conexpkPT_divorced(age,i,bignode,apnode,nodek)
      conexpp_divorced(age,i,bignode,apnode,nodek) = conexppPT_divorced(age,i,bignode,apnode,nodek)
      LFP_divorced(age,i,bignode,apnode,nodek) = 2
    else if (valuekNT_divorced(age,i,bignode,apnode,nodek) >= valuekFT_divorced(age,i,bignode,apnode,nodek) .and. &
             valuekNT_divorced(age,i,bignode,apnode,nodek) >= valuekPT_divorced(age,i,bignode,apnode,nodek)) then
      valuek_divorced(age,i,bignode,apnode,nodek) = valuekNT_divorced(age,i,bignode,apnode,nodek)
      valuep_divorced(age,i,bignode,apnode,nodek) = valuepNT_divorced(age,i,bignode,apnode,nodek)
      conexpk_divorced(age,i,bignode,apnode,nodek) = conexpkNT_divorced(age,i,bignode,apnode,nodek)
      conexpp_divorced(age,i,bignode,apnode,nodek) = conexppNT_divorced(age,i,bignode,apnode,nodek)
      LFP_divorced(age,i,bignode,apnode,nodek) = 1
    endif

  end do !end loop over insurance states

  !parent makes insurance decision in first period - only if healthy
  if (age==1) then 

    if (ni==3) then !public option
		if (valuep_divorced(age,1,bignode,apnode,nodek) .gt. valuep_divorced(age,2,bignode,apnode,nodek) &
			.and. valuep_divorced(age,1,bignode,apnode,nodek) .gt. valuep_divorced(age,3,bignode,apnode,nodek) &
			.and. health==1) then
		  insurance_divorced(bignode,apnode,nodek)=1
		elseif (valuep_divorced(age,3,bignode,apnode,nodek) .gt. valuep_divorced(age,1,bignode,apnode,nodek) &
			.and. valuep_divorced(age,3,bignode,apnode,nodek) .gt. valuep_divorced(age,2,bignode,apnode,nodek) &
			.and. health==1) then
		  insurance_divorced(bignode,apnode,nodek)=3 !public option
		else
		  insurance_divorced(bignode,apnode,nodek)=2		
		endif
	else 
		if (valuep_divorced(age,1,bignode,apnode,nodek) .gt. valuep_divorced(age,2,bignode,apnode,nodek) &
			.and. health==1) then
		  insurance_divorced(bignode,apnode,nodek)=1
		else
		  insurance_divorced(bignode,apnode,nodek)=2
		endif
	endif
	
  endif

 end subroutine solvevalfun_divorced
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine solvevalfunp_divorced(asschoicep,asschoicek,conk_div,lfp,ins,valp_div,conp_div,nodek)

  implicit none

  real(prec) asschoicep, asschoicek
  real(prec) conp_div, conk_div
  real(prec) valp_div
  integer lfp,ins
  integer medicaid

  integer nodek

  real(prec) cashp

  if (assetp(apnode,age)+incomep(wpnode,age)-hcost(ins,health,2)>cminp+mcdcash(health,2)) then
	cashp=assetp(apnode,age)+incomep(wpnode,age)-hcost(ins,health,2)
	medicaid=2
	if (ins==1) then
		cashp=cashp-premium
	elseif (ins==3) then
	    cashp=cashp-premiumPO
	endif
  else
	cashp=cminp+mcdcash(health,2)
	if (ins==1) then
		medicaid=1
		!medicaid=2  
		cashp=cashp-premium
	elseif (ins==3) then
		medicaid=1
		cashp=cashp-premiumPO
	else
		medicaid=1
	endif
  endif

  conp_div = cashp-asschoicep

  if (conp_div<0.0 .or. conp_div>cashp) then
    valp_div=-999999999d0
  else
    valp_div = valfun_p_divorced(asschoicep,asschoicek,conp_div,conk_div,lfp,ins,nodek,medicaid)
  endif

end subroutine solvevalfunp_divorced
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!solve KID value function in divorce
subroutine solvevalfunk_divorced(asschoicep,asschoicek,ins,valk_div,LFP_div,conk_div,nodek)
  
  implicit none

  real(prec) cashk
  real(prec) valk_div
  real(prec) leisure
  real(prec) conk_div
  integer LFP_div

  real(prec) asschoicep, asschoicek
  integer ins

  integer nodek
  integer wknode,aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  leisure=get_leisure(workhrs(LFP_div),famhrs(1))
  if (LFP_div==1) then 
    cashk=max(assetk(aknode,age) +nonwageinc(wknode)                              ,cmink)
  else if (LFP_div==2) then
    cashk=max(assetk(aknode,age)+PTwagefrac*incomek(wknode,age)+nonwageinc(wknode),cmink)
  else
    cashk=max(assetk(aknode,age)+           incomek(wknode,age)+nonwageinc(wknode),cmink)
  endif
  conk_div=cashk-asschoicek
  if (conk_div<0.0 .or. conk_div>cashk) then
    valk_div = -999999999d0
  else
    valk_div = valfun_k_divorced(asschoicep,asschoicek,conk_div,leisure,ins,nodek)
  endif

end subroutine solvevalfunk_divorced
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
subroutine solvevalfun_married(nodek)
  implicit none

  real(prec) cashp0, cashp100
  real(prec) cashkNT, cashkPT, cashkFT

  real(prec) cash,cashp
 
  real(prec) leisure

  real(prec) val
  real(prec) valtotNT0,valtotNT100,valtotPT0,valtotPT100,valtotFT0,valtotFT100

  real(prec) valtot, valp, valk
  real(prec) contot, conp, conk
  
  integer familycare
  integer insured
  integer medicaid0, medicaid100
  
  real(prec) ax, bx, cx


  integer k !counter for asset share (k) options

  integer nodek, wknode, aknode

  integer, parameter :: pastcmin=2 
  integer apnode_rip, aknode_rip
  integer nodep_rip, nodek_rip

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  !ax: minimum possible expenditure at time t
  ax = minass/2.d0

  if (apnode>pastcmin .and. aknode>pastcmin .and. aknode.ne.na) then

    !figure out the asset nodes to use
    if (apnode+aknode<=na+pastcmin) then
      apnode_rip=pastcmin
      aknode_rip=apnode+aknode-apnode_rip
    else
      aknode_rip=na
      apnode_rip=apnode+aknode-aknode_rip
    endif
    nodek_rip=(wknode-1)*na+aknode_rip

    !fill out all matrices
    do insured=1,ni
      valueNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valueNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaNT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaNT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valueNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valueNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaNT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaNT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)

      valuePT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuePT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepPT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepPT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekPT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekPT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppPT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppPT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkPT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkPT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaPT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaPT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuePT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuePT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepPT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepPT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekPT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekPT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppPT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppPT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkPT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkPT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaPT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaPT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)

      valueFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valueFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaFT0_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaFT0_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valueFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valueFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuepFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuepFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuekFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuekFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexppFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexppFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpkFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpkFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)         
      kappaFT100_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappaFT100_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)

      value_married(age,insured,bignode,mnode,apnode,nodek)= &
        value_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuep_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)= &
        valuek_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpp_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)= &
        conexpk_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)= &
        kappa_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)= &
        LFP_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
      famcare_married(age,insured,bignode,mnode,apnode,nodek)= &
        famcare_married(age,insured,bignode,mnode,apnode_rip,nodek_rip)
    enddo

    if (age==1) insurance_married(bignode,mnode,apnode,nodek)=insurance_married(bignode,mnode,apnode_rip,nodek_rip)

  !first time with this total amount of assets...
  else

  do insured=1,ni

	if (assetp(apnode,age)+incomep(wpnode,age)-hcost(insured,health,2)>cminp+mcdcash(health,2)) then
		cashp0=assetp(apnode,age)+incomep(wpnode,age)-hcost(insured,health,2) !0% family care
		medicaid0=2
		if (insured==1) then
			cashp0=cashp0-premium
		elseif (insured==3) then
			cashp0=cashp0-premiumPO
		endif
	else
		cashp0=cminp+mcdcash(health,2)
		if (insured==1) then
			!medicaid0=2  
			medicaid0=1  
			cashp0=cashp0-premium
		elseif (insured==3) then
			medicaid0=1  
			cashp0=cashp0-premiumPO			
		else
			medicaid0=1
		endif
	endif
	
	if (assetp(apnode,age)+incomep(wpnode,age)-hcost(insured,health,1)>cminp+mcdcash(health,1)) then
		cashp100=assetp(apnode,age)+incomep(wpnode,age)-hcost(insured,health,1) !100% family care
		medicaid100=2
		if (insured==1) then
			cashp100=cashp100-premium
		elseif (insured==3) then
			cashp100=cashp100-premiumPO
		endif
	else
		cashp100=cminp+mcdcash(health,1)
		medicaid100=1
		if (insured==1) then
			cashp100=cashp100-premium
		elseif (insured==3) then
			cashp100=cashp100-premiumPO			
		else
		endif
	endif  

    !compute values conditional on kid labor supply (0 and FT for now)

    !*******************************************
    !1. not working
      
      cashkNT=max(assetk(aknode,age)+0.0*incomek(wknode,age)+nonwageinc(wknode),cmink)
      
      !loop through care options (0% family, 100% family for now)
      
      !a. 0%
      familycare=0
      leisure=get_leisure(workhrs(1),famhrs(1))
      cash=cashp0+cashkNT
      cashp=cashp0
      cx=cash
      bx = cx*0.6d0
      if (age .eq. nt) then
        valtotNT0=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid0)
        valueNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotNT0
        valuepNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
        valuekNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
        conexppNT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
        conexpkNT0_married(age,insured,bignode,mnode,apnode,nodek)=conk
      else
        valtotNT0=9999999
        do k=1,nk
          val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid0)
          if (val<valtotNT0) then
            valtotNT0=val
            valueNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotNT0
            valuepNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
            valuekNT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
            conexppNT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
            conexpkNT0_married(age,insured,bignode,mnode,apnode,nodek)=conk            
            kappaNT0_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
          endif
        end do
      endif

      !b. 100% (only if unhealthy and not no-informal-care-simulation)
      if (health==1 .or. simnoinf==1) then
        valtotNT100=9999999
        valueNT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuepNT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuekNT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
      else
        familycare=1
        leisure=get_leisure(workhrs(1),famhrs(health)) !costs same as FTwork for now  
        cash=cashp100+cashkNT
        cashp=cashp100
        cx=cash
        bx = cx*0.6d0
        if (age .eq. nt) then
          valtotNT100=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																														medicaid100)
          valueNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotNT100
          valuepNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
          valuekNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
          conexppNT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
          conexpkNT100_married(age,insured,bignode,mnode,apnode,nodek)=conk
        else
          do k=1,nk
            valtotNT100=9999999
            val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid100)
            if (val<valtotNT100) then
              valtotNT100=val
              valueNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotNT100
              valuepNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
              valuekNT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
              conexppNT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
              conexpkNT100_married(age,insured,bignode,mnode,apnode,nodek)=conk            
              kappaNT100_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
            endif
          end do   
        endif  
      endif
    !*******************************************

    !*******************************************
    !1. part-time
      
      cashkPT=max(assetk(aknode,age)+PTwagefrac*incomek(wknode,age)+nonwageinc(wknode),cmink)
      
      !loop through care options (0% family, 100% family for now)
      
      !a. 0%
      familycare=0
      leisure=get_leisure(workhrs(2),famhrs(1))
      cash=cashp0+cashkPT
      cashp=cashp0
      cx=cash
      bx = cx*0.6d0
      if (age .eq. nt) then
        valtotPT0=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid0)
        valuePT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotPT0
        valuepPT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
        valuekPT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
        conexppPT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
        conexpkPT0_married(age,insured,bignode,mnode,apnode,nodek)=conk
      else
        valtotPT0=9999999
        do k=1,nk
          val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid0)
          if (val<valtotPT0) then
            valtotPT0=val
            valuePT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotPT0
            valuepPT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
            valuekPT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
            conexppPT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
            conexpkPT0_married(age,insured,bignode,mnode,apnode,nodek)=conk            
            kappaPT0_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
          endif
        end do
      endif

      !b. 100% (only if unhealthy and not no-informal-care-simulation)
      if (health==1 .or. simnoinf==1) then
        valtotPT100=9999999
        valuePT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuepPT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuekPT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
      else
        familycare=1
        leisure=get_leisure(workhrs(2),famhrs(health))   
        cash=cashp100+cashkPT
        cashp=cashp100
        cx=cash
        bx = cx*0.6d0
        if (age .eq. nt) then
          val=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																												medicaid100)
          valuePT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotPT100
          valuepPT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
          valuekPT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
          conexppPT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
          conexpkPT100_married(age,insured,bignode,mnode,apnode,nodek)=conk
        else
          do k=1,nk
            valtotPT100=9999999
            val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid100)
            if (val<valtotPT100) then
              valtotPT100=val
              valuePT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotPT100
              valuepPT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
              valuekPT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
              conexppPT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
              conexpkPT100_married(age,insured,bignode,mnode,apnode,nodek)=conk            
              kappaPT100_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
            endif
          end do   
        endif  
      endif
    !*******************************************    

!     !*******************************************
!     !2. full-time (only if necessary to compute)
      
      cashkFT=max(assetk(aknode,age)+1.0*incomek(wknode,age)+nonwageinc(wknode),cmink)
      
      !loop through care options (0% family, 100% family for now)
      
      !a. 0%
      familycare=0
      leisure=get_leisure(workhrs(3),famhrs(1))
      cash=cashp0+cashkFT
      cashp=cashp0
      cx=cash
      bx = cx*0.6d0
      if (age .eq. nt) then
        valtotFT0=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																														medicaid0)
        valueFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotFT0
        valuepFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
        valuekFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
        conexppFT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
        conexpkFT0_married(age,insured,bignode,mnode,apnode,nodek)=conk
      else
        valtotFT0=9999999
        do k=1,nk
          val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid0)
          if (val<valtotFT0) then
            valtotFT0=val
            valueFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valtotFT0
            valuepFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valp
            valuekFT0_married(age,insured,bignode,mnode,apnode,nodek)=-valk
            conexppFT0_married(age,insured,bignode,mnode,apnode,nodek)=conp
            conexpkFT0_married(age,insured,bignode,mnode,apnode,nodek)=conk            
            kappaFT0_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
          endif
        end do
      endif

      !b. 100% (only if unhealthy and not no-informal-care-simulation)
      if (health==1 .or. simnoinf==1) then
        valtotFT100=9999999
        valueFT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuepFT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
        valuekFT100_married(age,insured,bignode,mnode,apnode,nodek)=-99999999
      else
        familycare=1
        leisure=get_leisure(workhrs(3),famhrs(health)) 
        cash=cashp100+cashkFT
        cashp=cashp100
        cx=cash
        bx = cx*0.6d0
        if (age .eq. nt) then
          valtotFT100=golden_cons(ax,bx,cx,valfun_married,contot,0.5d0,leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																														medicaid100)
          valueFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotFT100
          valuepFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
          valuekFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
          conexppFT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
          conexpkFT100_married(age,insured,bignode,mnode,apnode,nodek)=conk
        else
          do k=1,nk
            valtotFT100=9999999
            val=golden_cons(ax,bx,cx,valfun_married,contot,kgrid(k),leisure,familycare,insured,conp,conk,valp,valk,cash,nodek, &
																													medicaid100)
            if (val<valtotFT100) then
              valtotFT100=val
              valueFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valtotFT100
              valuepFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valp
              valuekFT100_married(age,insured,bignode,mnode,apnode,nodek)=-valk
              conexppFT100_married(age,insured,bignode,mnode,apnode,nodek)=conp
              conexpkFT100_married(age,insured,bignode,mnode,apnode,nodek)=conk            
              kappaFT100_married(age,insured,bignode,mnode,apnode,nodek)=kgrid(k)
            endif
          end do   
        endif  
      endif
    !*******************************************

    valtot=min(valtotNT0,valtotNT100,valtotPT0,valtotPT100,valtotFT0,valtotFT100) 
    if (valtot==valtotNT0) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valueNT0_married(age,insured,bignode,mnode,apnode,nodek)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=valuepNT0_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=valuekNT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppNT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkNT0_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=kappaNT0_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=1
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=0
    else if (valtot==valtotNT100) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valueNT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuepNT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuekNT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppNT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkNT100_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=&
        kappaNT100_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=1
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=1
    else if (valtot==valtotPT0) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valuePT0_married(age,insured,bignode,mnode,apnode,nodek) 
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=valuepPT0_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=valuekPT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppPT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkPT0_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=&
        kappaPT0_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=2
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=0      
    else if (valtot==valtotPT100) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valuePT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuepPT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuekPT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppPT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkPT100_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=kappaPT100_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=2
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=1
    else if (valtot==valtotFT0) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valueFT0_married(age,insured,bignode,mnode,apnode,nodek)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=valuepFT0_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=valuekFT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppFT0_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkFT0_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=kappaFT0_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=3
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=0    
    else if (valtot==valtotFT100) then
      value_married(age,insured,bignode,mnode,apnode,nodek)=valueFT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuep_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuepFT100_married(age,insured,bignode,mnode,apnode,nodek)
      valuek_married(age,insured,bignode,mnode,apnode,nodek)=&
        valuekFT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpp_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexppFT100_married(age,insured,bignode,mnode,apnode,nodek)
      conexpk_married(age,insured,bignode,mnode,apnode,nodek)=&
        conexpkFT100_married(age,insured,bignode,mnode,apnode,nodek)
      kappa_married(age,insured,bignode,mnode,apnode,nodek)=kappaFT100_married(age,insured,bignode,mnode,apnode,nodek)
      LFP_married(age,insured,bignode,mnode,apnode,nodek)=3
      famcare_married(age,insured,bignode,mnode,apnode,nodek)=1
    else
      print*,'PROBLEM ASSIGNING MARRIAGE CONDITIONAL VALUE'
    endif

  enddo !end loop through insurance/no insurance

  !choose insurance or no insurance in first period
  if (age .eq. 1) then 
  
	if (ni==3) then
	    if (value_married(age,1,bignode,mnode,apnode,nodek)>=value_married(age,2,bignode,mnode,apnode,nodek) &
			.and. value_married(age,1,bignode,mnode,apnode,nodek)>=value_married(age,3,bignode,mnode,apnode,nodek) &
			.and. health==1) then
		  insurance_married(bignode,mnode,apnode,nodek)=1
	    elseif (value_married(age,3,bignode,mnode,apnode,nodek)>=value_married(age,1,bignode,mnode,apnode,nodek) &
				.and. value_married(age,3,bignode,mnode,apnode,nodek)>=value_married(age,2,bignode,mnode,apnode,nodek) &
				.and. health==1) then
		  insurance_married(bignode,mnode,apnode,nodek)=3
		else
		  insurance_married(bignode,mnode,apnode,nodek)=2
		endif
	else 
		if (value_married(age,1,bignode,mnode,apnode,nodek)>=value_married(age,2,bignode,mnode,apnode,nodek) &
			.and. health==1) then
		  insurance_married(bignode,mnode,apnode,nodek)=1
		else
		  insurance_married(bignode,mnode,apnode,nodek)=2
		endif
	endif
  endif

  endif !over whether calculations already done for that amount of total assets

end subroutine solvevalfun_married

!***************************************
real(prec) function valfun_married(conexp,kappaamt,leisure,famcare,insured, cp,ck,vp,vk,cash,nodek,medicaid)

  real(prec) conexp, cp, ck
  real(prec) kappaamt
  real(prec) leisure
  real(prec) vp, vk

  real(prec) cash

  integer famcare
  integer insured
  integer medicaid

  integer nodek

  ck=get_ck(conexp,leisure)
  cp=get_cp(conexp,ck)

  vk=valfun_k_married(conexp,ck,kappaamt,leisure,insured,cash,nodek)
  vp=valfun_p_married(conexp,cp,kappaamt,famcare,insured,leisure,cash,nodek,medicaid)
  valfun_married = mweight*vp + (1.0-mweight)*vk

end function valfun_married

real(prec) function valfun_p_married(conexp,cp,kappaamt,famcare,insured,leisure,cash,nodek,medicaid)

  real(prec) conexp, cp
  real(prec) leisure
  real(prec) kappaamt
  real(prec) utnow, utnowp, utnowk
  real(prec) saving_tot, saving_p, saving_k
  real(prec) node_interp_p, node_interp_k, node_interp_ktot
  real(prec) wgt1, wgt2
  real(prec) xpts(2),lb(2),ub(2)
  real(prec) valtom, valtombit
  real(prec) valkdead_tom, valkdead_tombit
  real(prec) cashtom,estim1,estim2

  real(prec) valterm
  real(prec) valpdead_tom

  real(prec) cash

  integer famcare
  integer insured
  integer medicaid
  integer mina_p, mina_k, mina_ktot
  integer minp
  integer bignodeto
  integer nodetok
  integer wwp, wwk
  integer hh
  integer llp, llk

  integer nodek, wknode, aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  wgt1=0.d0
  wgt2=0.d0

  utnowk = utility_married(conexp-cp,leisure,0,2,2) !consumption, leisure, famcare (placeholder), kid, medicaid (placeholder)
  utnowp = utility_married(cp,0.0d0,famcare,1,medicaid) !consumption, leisure (placeholder), famcare, parent, medicaid
  utnow  = utnowp + altruism*utnowk

  !find next period asset stock
  saving_tot = r*(cash - conexp)
  saving_p=max(minass,kappaamt*saving_tot)
  saving_k=max(minass,(1.0-kappaamt)*saving_tot)

  if (ageto.ge.nt+1) then !if last period (nt)
      
    valterm = utilitykterminal(saving_tot,wknode) !kid terminal value function (parent cares about this through altruism - like bequest motive)
    valterm = valterm*altruism
    valfun_p_married = - (utnow + (disc**(-1.d0))*valterm)

  else     !if not last period

    !parent savings:
    call find_ass(1, saving_p, ageto, mina_p, wgt1,wgt2)
    node_interp_p=wgt1*real(mina_p)+wgt2*real(mina_p+1)
    !kid savings:
    call find_ass(2, saving_k, ageto, mina_k, wgt1,wgt2)
    node_interp_k=wgt1*real(mina_k)+wgt2*real(mina_k+1)
    !total savings (relevant if parent dies):
    call find_ass(2, saving_tot, ageto, mina_ktot,wgt1,wgt2)
    node_interp_ktot=wgt1*real(mina_ktot)+wgt2*real(mina_ktot+1)

    !define parent's asset node
    lb(1)=mina_p
    ub(1)=mina_p+1
    xpts(1)=node_interp_p

    valtom = 0.d0
    valkdead_tom = 0.d0
    !loop through kid future income states (not parent cuz parent income is deterministic)
    do wwk = maxwn(wknode,ageto)+1 ,minwn(wknode,ageto)-1
      
      !define kid's income*asset (compact) node
      nodetok = (wwk-1)*na
      lb(2)=nodetok+mina_k
      ub(2)=nodetok+mina_k+1                
      xpts(2)=nodetok+node_interp_k

      valkdead_tombit=interp_deadparent(ageto,nodetok,mina_ktot,node_interp_ktot,valuek_deadparent)
      valkdead_tombit=valkdead_tombit*wprobk(wwk,wknode,ageto)
      valkdead_tom=valkdead_tom+valkdead_tombit

      do hh=1,nh

        bignodeto = (hh-1)*nwp + wpnode

          !bilinear interpolation: 2 continuous variables are (1) parent cash, (2) kid cash
          valtombit=bilinear(xpts,lb,ub, &
                      valuep(ageto,insured,bignodeto,mnode,int(lb(1)):int(ub(1)),int(lb(2)):int(ub(2))))

          !weight by prob of ww (income) at t+1 given wnode at t (just for k)
          valtombit=valtombit*wprobk(wwk,wknode,ageto)
          !weight by prob of hh (health) at t+1 given hnode at t (just for p)
          valtombit=valtombit*hprob(hh,health,wpnode,ageto)

          valtom = valtom + valtombit

      enddo !end loop over future health states
    enddo   !end loop over future kid income states

    valpdead_tom=altruism*valkdead_tom

    valfun_p_married = - (utnow + (disc**(-1.d0))*(       survprob(health,wpnode,ageto)*valtom + &
                                                   (1.0d0-survprob(health,wpnode,ageto))*valpdead_tom))

  endif

end function valfun_p_married


real(prec) function valfun_k_married(conexp,ck,kappaamt,leisure,insured,cash,nodek)

  real(prec) conexp, ck
  real(prec) kappaamt
  real(prec) leisure
  real(prec) utnow
  real(prec) vterm
  real(prec) saving_tot, saving_p, saving_k
  real(prec) node_interp_p, node_interp_k, node_interp_ktot
  real(prec) pwgt_interp
  real(prec) wgt1, wgt2
  real(prec) xpts(2),lb(2),ub(2)
  real(prec) valtom, valtombit
  real(prec) cashtom,estim1,estim2
  real(prec) valterm

  real(prec) cash

  real(prec) valkdead_tom, valkdead_tombit

  integer insured
  integer mina_p, mina_k, mina_ktot
  integer minp
  integer nodetop, nodetok
  integer wwp, wwk
  integer hh
  integer llp, llk
  integer bignodeto

  integer nodek, wknode, aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  wgt1=0.d0
  wgt2=0.d0

  utnow = utility_married(ck,leisure,0,2,2) !consumption, leisure, famcare (placeholder),kid, medicaid (placeholder)

  !find savings (for next period or for kid terminal condition)
  saving_tot = r*(cash - conexp)
  saving_p=max(minass,kappaamt*saving_tot)
  saving_k=max(minass,(1.0-kappaamt)*saving_tot)

  if (ageto.ge.nt+1) then !if last period (nt)

    valterm=utilitykterminal(saving_tot,wknode)
    valfun_k_married = -(utnow + (disc**(-1.d0))*valterm)

  else     !if not last period

    !search optimal savings division kappa (fraction going to parent)

    !parent savings:
    call find_ass(1, saving_p, ageto, mina_p, wgt1,wgt2)
    node_interp_p=wgt1*real(mina_p)+wgt2*real(mina_p+1)
    !kid savings:
    call find_ass(2, saving_k, ageto, mina_k, wgt1,wgt2)
    node_interp_k=wgt1*real(mina_k)+wgt2*real(mina_k+1)
    !total savings (relevant if parent dies):
    call find_ass(2, saving_tot, ageto, mina_ktot,wgt1,wgt2)
    node_interp_ktot=wgt1*real(mina_ktot)+wgt2*real(mina_ktot+1)

    !define parent's asset node
    lb(1)=mina_p
    ub(1)=mina_p+1
    xpts(1)=node_interp_p
        
    !loop through kid future income states (not parent cuz parent income is deterministic)
    valtom = 0.d0
    valkdead_tom = 0.d0
    do wwk = maxwn(wknode,ageto)+1 ,minwn(wknode,ageto)-1
      
      !define income*asset (compact) node
      nodetok = (wwk-1)*na
      lb(2)=nodetok+mina_k
      ub(2)=nodetok+mina_k+1                
      xpts(2)=nodetok+node_interp_k

      valkdead_tombit=interp_deadparent(ageto,nodetok,mina_ktot,node_interp_ktot,valuek_deadparent)
      valkdead_tombit=valkdead_tombit*wprobk(wwk,wknode,ageto)
      valkdead_tom=valkdead_tom+valkdead_tombit

      do hh=1,nh

        bignodeto = (hh-1)*nwp + wpnode

          !trilinear interpolation: 3 continuous variables are (1) parent cash, (2) kid cash, (3) pareto weight
          valtombit=bilinear(xpts,lb,ub, &
                      valuek(ageto,insured,bignodeto,mnode,int(lb(1)):int(ub(1)),int(lb(2)):int(ub(2))))

          !weight by prob of ww (income) at t+1 given wnode at t (just for k)
          valtombit=valtombit*wprobk(wwk,wknode,ageto)
          !weight by prob of hh (health) at t+1 given hnode at t (just for p)
          valtombit=valtombit*hprob(hh,health,wpnode,ageto)

          valtom = valtom + valtombit

      enddo !end loop over future health states
    enddo   !end loop over future kid income states

    valfun_k_married = - (utnow + (disc**(-1.d0))*(       survprob(health,wpnode,ageto)*valtom + &
                                                   (1.0d0-survprob(health,wpnode,ageto))*valkdead_tom)) 

  endif

end function valfun_k_married
!***************************************

!***************************************
real(prec) function valfun_p_divorced(asschoicep,asschoicek,conexp_p,conexp_k,LFP,insured,nodek,medicaid)

  real(prec) asschoicep,asschoicek
  real(prec) conexp_p,conexp_k
  integer LFP
  integer insured
  integer medicaid

  real(prec) leisurek
  real(prec) saving_p, saving_k, saving_tot
  real(prec) utnowp,utnowk,utnow
  real(prec) node_interp_p, node_interp_k, node_interp_ktot
  integer mina_p, mina_k, mina_ktot
  integer nodetop, nodetok
  real(prec) xpts(2),lb(2),ub(2)

  real(prec) valtom, valtombit
  real(prec) valkdead_tom, valkdead_tombit
  real(prec) valpdead_tom

  real(prec) wgt1,wgt2
  real(prec) valterm
  integer wwk,hh
  integer bignodeto

  integer nodek, wknode, aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  wgt1=0.d0
  wgt2=0.d0

  !calculate flow utility this period
  leisurek=get_leisure(workhrs(LFP),0.d0)
  utnowk = utility_divorced(conexp_k,leisurek,2,2) !consumption, leisure, kid, not on medicaid (placeholder)
  utnowp = utility_divorced(conexp_p,0.0d0   ,1,medicaid) !consumption, leisure (placeholder), parent, whether on medicaid
  utnow  = utnowp + altruism*utnowk

  !calculate kid terminal value
  saving_k = r*asschoicek
  if (saving_k.le.minass) saving_k = minass
  saving_p = r*asschoicep
  if (saving_p.le.minass) saving_p = minass
  saving_tot = r*(asschoicek+asschoicep)


  !if last period (nt)
  if (ageto.ge.nt+1) then

    valterm=utilitykterminal(saving_k+saving_p,wknode)    
    valterm=altruism*valterm

    valfun_p_divorced = utnow + valterm
  
  !if not last period
  else

    !if not in final period, find next period asset stock on grid & calculate expected value (given alive next period)

    !parent savings:
    call find_ass(1, saving_p, ageto, mina_p, wgt1,wgt2)
    node_interp_p=wgt1*real(mina_p)+wgt2*real(mina_p+1)
    !kid savings:
    call find_ass(2, saving_k, ageto, mina_k, wgt1,wgt2)
    node_interp_k=wgt1*real(mina_k)+wgt2*real(mina_k+1)
    !total savings (relevant if parent dies):
    call find_ass(2, saving_tot, ageto, mina_ktot,wgt1,wgt2)
    node_interp_ktot=wgt1*real(mina_ktot)+wgt2*real(mina_ktot+1)

    !define parent's asset node
    lb(1)=mina_p
    ub(1)=mina_p+1
    xpts(1)=node_interp_p

    valtom = 0.d0
    valkdead_tom = 0.d0
    !loop through kid future income states (not parent cuz parent income is deterministic)
    do wwk = maxwn(wknode,ageto)+1 ,minwn(wknode,ageto)-1

      !define kid's income*asset (compact) node
      nodetok = (wwk-1)*na
      lb(2)=nodetok+mina_k
      ub(2)=nodetok+mina_k+1                
      xpts(2)=nodetok+node_interp_k


      valkdead_tombit=interp_deadparent(ageto,nodetok,mina_ktot,node_interp_ktot,valuek_deadparent)
      valkdead_tombit=valkdead_tombit*wprobk(wwk,wknode,ageto)
      valkdead_tom=valkdead_tom+valkdead_tombit

      do hh=1,nh

          bignodeto = (hh-1)*nwp + wpnode

          !bilinear interpolation: 2 continuous variables are (1) parent cash, (2) kid cash
          valtombit=bilinear(xpts,lb,ub, &
                      valuep_divorced(ageto,insured,bignodeto,int(lb(1)):int(ub(1)),int(lb(2)):int(ub(2))))

          !weight by prob of ww (income) at t+1 given wnode at t (just for k)
          valtombit=valtombit*wprobk(wwk,wknode,ageto)
          !weight by prob of hh (health) at t+1 given hnode at t (just for p)
          valtombit=valtombit*hprob(hh,health,wpnode,ageto)

          valtom = valtom + valtombit

      enddo !end loop over future health states
    enddo   !end loop over future kid income states

    valpdead_tom=altruism*valkdead_tom

    valfun_p_divorced = utnow + (disc**(-1.d0))*(        survprob(health,wpnode,ageto)*valtom + &
                                                 (1.0d0-survprob(health,wpnode,ageto))*valpdead_tom)

  endif

end function valfun_p_divorced
!***************************************


!***************************************
real(prec) function valfun_k_divorced(asschoicep,asschoicek,conexp,leisure,insured,nodek)

  real(prec) asschoicep,asschoicek
  real(prec) conexp
  real(prec) leisure
  integer insured

  real(prec) utnow
  real(prec) saving_p, saving_k, saving_tot
  real(prec) node_interp_p, node_interp_k, node_interp_ktot
  integer mina_p, mina_k, mina_ktot
  integer nodetop, nodetok
  real(prec) xpts(2),lb(2),ub(2)

  real(prec) valtom, valtombit

  real(prec) valkdead_tom, valkdead_tombit
  real(prec) valpdead_tom

  real(prec) wgt1,wgt2
  real(prec) valterm
  integer wwk,hh
  integer bignodeto

  integer nodek, wknode, aknode

  wknode=ceiling(real(nodek)/real(na))
  aknode=nodek-(wknode-1)*na 

  wgt1=0.d0
  wgt2=0.d0

  utnow =   utility_divorced(conexp,leisure,2,2) !consumption, leisure, kid, not on medicaid (placeholder)

  !calculate kid terminal value
  saving_k = r*asschoicek
  if (saving_k.le.minass) saving_k = minass
  saving_p = r*asschoicep
  if (saving_p.le.minass) saving_p = minass
  saving_tot = r*(asschoicek+asschoicep)

  !if last period (nt)
  if (ageto.ge.nt+1) then

    valterm=utilitykterminal(saving_tot,wknode)     

    valfun_k_divorced = utnow + (disc**(-1.d0))*valterm
  
  !if not last period
  else

    !if not in final period, find next period asset stock on grid & calculate expected value (given alive next period)

    !parent savings:
    call find_ass(1, saving_p, ageto, mina_p, wgt1,wgt2)
    node_interp_p=wgt1*real(mina_p)+wgt2*real(mina_p+1)
    !kid savings:
    call find_ass(2, saving_k, ageto, mina_k, wgt1,wgt2)
    node_interp_k=wgt1*real(mina_k)+wgt2*real(mina_k+1)
    !total savings (relevant if parent dies):
    call find_ass(2, saving_tot, ageto, mina_ktot,wgt1,wgt2)
    node_interp_ktot=wgt1*real(mina_ktot)+wgt2*real(mina_ktot+1)

    !define parent's asset node
    lb(1)=mina_p
    ub(1)=mina_p+1
    xpts(1)=node_interp_p

    valtom = 0.d0
    valkdead_tom = 0.d0
    !loop through kid future income states (not parent cuz parent income is deterministic)
    do wwk = maxwn(wknode,ageto)+1 ,minwn(wknode,ageto)-1
      !define kid's income*asset (compact) node
      nodetok = (wwk-1)*na
      lb(2)=nodetok+mina_k
      ub(2)=nodetok+mina_k+1                
      xpts(2)=nodetok+node_interp_k

      valkdead_tombit=interp_deadparent(ageto,nodetok,mina_ktot,node_interp_ktot,valuek_deadparent)
      valkdead_tombit=valkdead_tombit*wprobk(wwk,wknode,ageto)
      valkdead_tom=valkdead_tom+valkdead_tombit

      do hh=1,nh

          bignodeto = (hh-1)*nwp + wpnode

          !bilinear interpolation: 2 continuous variables are (1) parent cash, (2) kid cash
          valtombit=bilinear(xpts,lb,ub, &
                      valuek_divorced(ageto,insured,bignodeto,int(lb(1)):int(ub(1)),int(lb(2)):int(ub(2))))

          !weight by prob of ww (income) at t+1 given wnode at t (just for k)
          valtombit=valtombit*wprobk(wwk,wknode,ageto)
          !weight by prob of hh (health) at t+1 given hnode at t (just for p)
          valtombit=valtombit*hprob(hh,health,wpnode,ageto)

          valtom = valtom + valtombit

      enddo !end loop over future health states
    enddo   !end loop over future kid income states

    valfun_k_divorced = utnow + (disc**(-1.d0))*(        survprob(health,wpnode,ageto)*valtom + &
                                                 (1.0d0-survprob(health,wpnode,ageto))*valkdead_tom)

  endif

end function valfun_k_divorced
!***************************************


end module valuefunction

